* This MC code generates data with partial observability and assesses the impact of revealing a percentage of 0 cases


clear
set more off

* This is homoskedastic PI spec model, no extra info
capture program drop partobs1_lf
program define partobs1_lf
        version 8.0
	args 	lnf uV2 uV3 uI3	
	tempvar pV pI 
	qui gen double `pI' = norm((`uI3')/sqrt(2))
	qui gen double `pV' = norm((`pI'*`uV3' + (1-`pI')*`uV2')/sqrt(2))

	qui replace `lnf' = ln(`pV') + ln(`pI')   if $ML_y1 ==1 
	qui replace `lnf' = ln(1-`pV'*`pI') if $ML_y1 == 0 | $ML_y1 == 2 | $ML_y1 == 3
end



* This is homoskedastic PI spec model, with some observ revealed
capture program drop partobs2_lf
program define partobs2_lf
        version 8.0
	args 	lnf uV2 uV3 uI3	
	tempvar pV pI 
	qui gen double `pI' = norm((`uI3')/sqrt(2))
	qui gen double `pV' = norm((`pI'*`uV3' + (1-`pI')*`uV2')/sqrt(2))

	qui replace `lnf' = ln(`pV') + ln(`pI')   if $ML_y1 == 1
	qui replace `lnf' = ln(1-`pV'*`pI') if $ML_y1 == 0
	qui replace `lnf' = ln(`pV') + ln(1-`pI') if $ML_y1 == 2
	qui replace `lnf' = ln(1-`pV') if $ML_y1 == 3	
end


capture program drop PartMC
program PartMC, rclass
	version 8.0
	args obser
	qui set obs `obser'

local seed = "`c(seed)'"

** DATA GENERATING PROCESS
* generate regresssors
	qui gen double x11 = rnormal()
	qui gen double x12 = rnormal()
	qui gen double x13 = rnormal()
	qui gen double x13c = uniform()
	qui gen double x23 = rnormal()
	qui gen double x23c = uniform()
	qui gen double xcom = uniform()

* generate errors for player 2
qui gen double e22 = rnormal(0, 1)
qui gen double e23 = rnormal(0, 1)

* 2's probability
qui gen double p2_res = normal((.25 + x23 -x23c + .5*xcom)/sqrt(2))


* generate errors for player 1
qui gen double e11 = rnormal(0, 1)
qui gen double e12 = rnormal(0, 1)

* 1's probability
qui gen double p1_att = normal((p2_res*(.5+1.5*x13-x13c -.5*xcom)+(1-p2_res)*(.3-.4*x12))/sqrt(2))

* B's DV, if fully observed
qui gen y2 = 0
qui replace y2 = 1 if  .25 + x23 - x23c + .5*xcom + e23 > e22

qui sum y2
scalar y2freq = r(mean)

* A's DV, if fully observed
qui gen y1 = 0
qui replace y1 = 1 if p2_res*(.5+1.5*x13-x13c -.5*xcom)+(1-p2_res)*(.3 -.4*x12) +e12 > e11 
qui sum y1
scalar y1freq = r(mean)


* fully observed, 3 outcomes coded
gen start2 = .
replace start2 = 1 if y1==1 & y2==1
replace start2 = 2 if y1==1 & y2==0
replace start2 = 3 if y1==0

* partially observed, only 1 and 0 outcomes
gen starta=.
replace starta = 1 if y1==1 & y2==1
replace starta = 0 if y1!=1 | y2!=1


* how many 0s exist originally
count if starta==0
scalar Arep0 = r(N)

count if starta==1
scalar Arep1 = r(N)

* random revelation
gen rand = uniform()
g rand2 = round(uniform(), .01) in 1
replace rand2 = rand2[1]

scalar rev = rand2[1]

gen start = starta
replace start=start2 if rand<=rand2


* how many 2's are revealed
count if rand<=rand2 & start2==2
scalar Srep2 = r(N)

* how many 3's are revealed
count if rand<=rand2 & start2==3
scalar Srep3 = r(N)

* how many 0s remain after 5% revelation
count if start==0
scalar Srep0 = r(N)


qui gen double prwar = p1_att*p2_res


* estimating the partial obs model, without any additional info

ml model lf partobs1_lf  (uV2:start = x12) (uV3:start = x13 x13c xcom) (uI3: start = x23 x23c xcom)
capture {
ml maximize, difficult iterate(50)
	matrix Hbprob = e(b)
	matrix HV = vecdiag(e(V))	
scalar converH = e(converged)
predict double u_12, eq(#1)
predict double u_13, eq(#2)
predict double u_23, eq(#3)
	qui gen double p2resp = normal((u_23)/sqrt(2))
	qui gen double p1attp = normal((p2resp*u_13+(1-p2resp)*u_12)/sqrt(2))
qui gen double prwarH = p1attp*p2resp
drop u_*  p1attp
}

if _rc != 0 {
	matrix Hbprob = (.,.,.,.,.,.,.,.,.,.)
	matrix HV = (.,.,.,.,.,.,.,.,.,.)
}
if _rc == 0 {
if converH==0{
	matrix Hbprob = (.,.,.,.,.,.,.,.,.,.)
	matrix HV = (.,.,.,.,.,.,.,.,.,.)
}
}


* estimating the partial obs model, with x% revelation
ml model lf partobs2_lf  (uV2:start = x12) (uV3:start = x13 x13c xcom) (uI3: start = x23 x23c xcom)
capture {
ml maximize, difficult iterate(50)
	matrix Hbprob2 = e(b)
	matrix HV2 = vecdiag(e(V))	
scalar converH2 = e(converged)
predict double u_12a, eq(#1)
predict double u_13a, eq(#2)
predict double u_23a, eq(#3)
	qui gen double p2resp2 = normal((u_23a)/sqrt(2))
	qui gen double p1attp2 = normal((p2resp2*u_13a+(1-p2resp2)*u_12a)/sqrt(2))
qui gen double prwarH2 = p1attp2*p2resp2
drop u_*  p1attp2
}

if _rc != 0 {
	matrix Hbprob2 = (.,.,.,.,.,.,.,.,.,.)
	matrix HV2 = (.,.,.,.,.,.,.,.,.,.)
}
if _rc == 0 {
if converH2==0{
	matrix Hbprob2 = (.,.,.,.,.,.,.,.,.,.)
	matrix HV2 = (.,.,.,.,.,.,.,.,.,.)
}
}



* estimating the partial obs model, full info
ml model lf partobs2_lf  (uV2:start2 = x12) (uV3:start2 = x13 x13c xcom) (uI3: start2 = x23 x23c xcom)
capture {
ml maximize, difficult iterate(50)
	matrix Hbprob5 = e(b)
	matrix HV5 = vecdiag(e(V))	
scalar converH5 = e(converged)
predict double u_12b, eq(#1)
predict double u_13b, eq(#2)
predict double u_23b, eq(#3)
	qui gen double p2resp3 = normal((u_23b)/sqrt(2))
	qui gen double p1attp3 = normal((p2resp3*u_13b+(1-p2resp3)*u_12b)/sqrt(2))
qui gen double prwarH5 = p1attp3*p2resp3
drop u_*  p1attp3
}

if _rc != 0 {
	matrix Hbprob5 = (.,.,.,.,.,.,.,.,.,.)
	matrix HV5 = (.,.,.,.,.,.,.,.,.,.)
}
if _rc == 0 {
if converH5==0{
	matrix Hbprob5 = (.,.,.,.,.,.,.,.,.,.)
	matrix HV5 = (.,.,.,.,.,.,.,.,.,.)
}
}





capture {
gen double rmseHom =  (p2_res- p2resp)^2
sum rmseHom
scalar diffHom = sqrt(r(mean))
scalar diffmaxHom = sqrt(r(max))
gen temperc = sqrt(rmseHom)/p2_res
centile temperc
scalar diffpercHom = r(c_1)
}
if _rc != 0{
scalar diffHom =.
scalar diffmaxHom =.
scalar diffpercHom = .
}
if Hbprob[1,1]==.{
scalar diffHom =.
scalar diffmaxHom =.
scalar diffpercHom = .
}


capture {
gen double rmseHet =  (p2_res- p2resp2)^2
sum rmseHet
scalar diffHet = sqrt(r(mean))
scalar diffmaxHet = sqrt(r(max))
gen temperc2 = sqrt(rmseHet)/p2_res
centile temperc2
scalar diffpercHet = r(c_1)
}
if _rc != 0{
scalar diffHet =.
scalar diffmaxHet =.
scalar diffpercHet = .
}
if Hbprob2[1,1]==.{
scalar diffHet =.
scalar diffmaxHet =.
scalar diffpercHet = .
}


capture {
gen double rmsePIBB =  (p2_res- p2resp3)^2
sum rmsePIBB
scalar diffPIBB = sqrt(r(mean))
scalar diffmaxPIBB = sqrt(r(max))
gen temperc5 = sqrt(rmsePIBB)/p2_res
centile temperc5
scalar diffpercPIBB = r(c_1)
}
if _rc != 0{
scalar diffPIBB =.
scalar diffmaxPIBB =.
scalar diffpercPIBB = .
}
if Hbprob5[1,1]==.{
scalar diffPIBB =.
scalar diffmaxPIBB =.
scalar diffpercPIBB = .
}




drop _all

	ret scalar rev = rev
	ret scalar y1freq = y1freq
	ret scalar y2freq = y2freq

	ret scalar diffHom = diffHom
	ret scalar diffmaxHom = diffmaxHom
	ret scalar diffpercHom = diffpercHom

	ret scalar diffHet = diffHet
	ret scalar diffmaxHet = diffmaxHet
	ret scalar diffpercHet = diffpercHet

	ret scalar diffPIBB = diffPIBB
	ret scalar diffmaxPIBB = diffmaxPIBB
	ret scalar diffpercPIBB = diffpercPIBB

	ret scalar Hbu12  = Hbprob[1,1]	
	ret scalar Hbu120  = Hbprob[1,2]	
	ret scalar Hbu13  = Hbprob[1,3]	
	ret scalar Hbu13c  = Hbprob[1,4]	
	ret scalar Hbu13cc  = Hbprob[1,5]
	ret scalar Hbu130  = Hbprob[1,6]	
	ret scalar Hbu23  = Hbprob[1,7]	
	ret scalar Hbu23c  = Hbprob[1,8]	
	ret scalar Hbu23cc  = Hbprob[1,9]
	ret scalar Hbu230  = Hbprob[1,10]	

	ret scalar Sbu12  = Hbprob2[1,1]	
	ret scalar Sbu120  = Hbprob2[1,2]	
	ret scalar Sbu13  = Hbprob2[1,3]	
	ret scalar Sbu13c  = Hbprob2[1,4]	
	ret scalar Sbu13cc  = Hbprob2[1,5]
	ret scalar Sbu130  = Hbprob2[1,6]	
	ret scalar Sbu23  = Hbprob2[1,7]	
	ret scalar Sbu23c  = Hbprob2[1,8]	
	ret scalar Sbu23cc  = Hbprob2[1,9]
	ret scalar Sbu230  = Hbprob2[1,10]	

	ret scalar BBBbu12  = Hbprob5[1,1]	
	ret scalar BBBbu120  = Hbprob5[1,2]	
	ret scalar BBBbu13  = Hbprob5[1,3]	
	ret scalar BBBbu13c  = Hbprob5[1,4]	
	ret scalar BBBbu13cc  = Hbprob5[1,5]
	ret scalar BBBbu130  = Hbprob5[1,6]	
	ret scalar BBBbu23  = Hbprob5[1,7]	
	ret scalar BBBbu23c  = Hbprob5[1,8]	
	ret scalar BBBbu23cc  = Hbprob5[1,9]
	ret scalar BBBbu230  = Hbprob5[1,10]	
	

	ret scalar Arep1 = Arep1

	
	ret scalar Arep2 = 0
	ret scalar Srep2 = Srep2
	
	ret scalar Arep3 = 0
	ret scalar Srep3 = Srep3
	
	ret scalar Arep0 = Arep0
	ret scalar Srep0 = Srep0
		
end

#delimit;
set seed 654321;
set more off;
simulate "PartMC 2000" 

	Reveal = r(rev)
	y1freq = r(y1freq)
	y2freq = r(y2freq)

	rmspH = r(diffHom)
	rmspS = r(diffHet)
	rmspBBB = r(diffPIBB)
	
	rmxspHmax = r(diffmaxHom)
	rmxspSmax = r(diffmaxHet)
	rmxspBBBmax = r(diffmaxPIBB)

	rmspHperc = r(diffpercHom)
	rmspSperc = r(diffpercHet)
	rmspBBBperc = r(diffpercPIBB)

	Arep1 = r(Arep1)

	Arep2 = r(Arep2)
	Srep2 = r(Srep2)
	
	Arep3 = r(Arep3)
	Srep3 = r(Srep3)
	
	Arep0 = r(Arep0)
	Srep0 = r(Srep0)
	
	
    	Hbu12  = r(Hbu12)
    	Sbu12  = r(Sbu12)
    	BBBbu12  = r(BBBbu12)
    	
		Hbu120  = r(Hbu120)
    	Sbu120  = r(Sbu120)
    	BBBbu120  = r(BBBbu120)

    	Hbu13  = r(Hbu13)
    	Sbu13  = r(Sbu13)
    	BBBbu13  = r(BBBbu13)
   	
		Hbu13c  = r(Hbu13c)
		Sbu13c  = r(Sbu13c)
		BBBbu13c  = r(BBBbu13c)
    	
		Hbu13cc  = r(Hbu13cc)
    	Sbu13cc  = r(Sbu13cc)
    	BBBbu13cc  = r(BBBbu13cc)
    	
		Hbu130  = r(Hbu130)
    	Sbu130  = r(Sbu130)
    	BBBbu130  = r(BBBbu130)

    	Hbu23  = r(Hbu23)
    	Sbu23  = r(Sbu23)
    	BBBbu23  = r(BBBbu23)
   	
		Hbu23c  = r(Hbu23c)
		Sbu23c  = r(Sbu23c)
		BBBbu23c  = r(BBBbu23c)
    	
		Hbu23cc  = r(Hbu23cc)
    	Sbu23cc  = r(Sbu23cc)
    	BBBbu23cc  = r(BBBbu23cc)
    	
		Hbu230  = r(Hbu230)
    	Sbu230  = r(Sbu230)
    	BBBbu230  = r(BBBbu230)

,reps(10000) dots saving(C:\Users\MB\Desktop\BasStoneISQ\MonteCarlo_N2000_I10000.dta) every(20) double replace;


* if you need to plot the figure 2 only, load MonteCarlo_N2000_I10000.dta and run below:

gen Reveal1 = (Arep0-Srep0)/Arep0

sum rmspH
gen rmspHplot = r(mean)
sum rmspBBB
gen rmspBBBplot = r(mean)
replace Reveal = round(Reveal, .01)
bysort Reveal: egen mnrmspS = mean(rmspS)

twoway (lpoly rmspS Reveal1, lcolor(black) degree(4)) (scatter mn* Reveal, mcolor(gray)), leg(off) ytitle("RMSE of Democracy") xtitle("Proportion of 0 Cases Revealed")